using System;
using System.Collections.Generic;
using dnAnalytics.LinearAlgebra;
using NUnit.Framework;

namespace dnAnalytics.Tests.LinearAlgebra
{
    [TestFixture]
    public class DiagonalMatrixTests
    {
        protected static readonly double[] diag3 = new double[] { -1.1, -2.2, -3.3 };
        protected static readonly double[] singular = new double[] { -3.3, -4.4, 0.0 };
        protected static readonly double[] diag5 = new double[] { -1.1, -2.2, -3.3, -4.4, -5.5 };
        protected static readonly double[] diag7 = new double[] { -1.1, -2.2, -3.3, -4.4, -5.5, -6.6, -7.7 };
        protected static readonly double[][] matrices = new double[][] { diag3, singular, diag5, diag7 };

        public DiagonalMatrix GetMatrix(double[] diag)
        {
            return new DiagonalMatrix(diag);
        }

        public DiagonalMatrix GetMatrix(int row, int col)
        {
            return new DiagonalMatrix(row, col);
        }

        public Matrix GetNonDiagonalMatrix(int row, int col)
        {
            return new DenseMatrix(row, col, 5.0);
        }

        public DenseVector GetVector(double[] vals)
        {
            return new DenseVector(vals);
        }

        public DenseVector GetVector(int size)
        {
            return new DenseVector(size);
        }

        [Test]
        public void EqualsTest()
        {
            Matrix A = GetMatrix(diag3);
            Matrix B = GetMatrix(diag3);
            Matrix C = A;
            Matrix D = GetMatrix(diag5);
            Matrix E = GetMatrix(singular);

            Assert.IsTrue(A.Equals(A));
            Assert.IsTrue(A.Equals(B));
            Assert.IsTrue(A.Equals(C));
            Assert.IsTrue(A.Equals(B) == B.Equals(A));
            Assert.IsTrue(A.Equals(C) == C.Equals(A));
            Assert.IsFalse(A.Equals(null));
            Assert.IsFalse(A.Equals(E));
            Assert.IsFalse(D.Equals(A));

            Assert.IsTrue(A == C);
            Assert.IsFalse(A != C);
            Assert.IsFalse(A == B);
            Assert.IsTrue(((IEquatable<Matrix>)A).Equals(B));
            Assert.IsTrue(((IEquatable<Matrix>)B).Equals(C));
        }

        [Test]
        public void ConstructorExceptions()
        {
            Assert.Throws<ArgumentException>(() => GetMatrix(0, 5));
            Assert.Throws<ArgumentException>(() => GetMatrix(5, 0));
        }

        [Test]
        public void ConstructorValue()
        {
            const int dim = 4;
            DiagonalMatrix matrix = new DiagonalMatrix(dim, 2.0);
            for (int row = 0; row < dim; row++)
            {
                for (int col = 0; col < dim; col++)
                {
                    if (row == col)
                    {
                        Assert.AreEqual(2.0, matrix[row, col]);
                    }
                    else
                    {
                        Assert.AreEqual(matrix[row, col], 0.0);
                    }
                }
            }
        }

        [Test]
        public void ConstructorDenseVector()
        {
            const int dim = 4;
            DenseVector dense = new DenseVector(dim);
            dense[0] = 2.0;
            dense[3] = 3.0;
            DiagonalMatrix matrix = new DiagonalMatrix(dense);
            for (int row = 0; row < dim; row++)
            {
                for (int col = 0; col < dim; col++)
                {
                    if (row == col)
                    {
                        Assert.AreEqual(dense[row], matrix[row, col]);
                    }
                    else
                    {
                        Assert.AreEqual(matrix[row, col], 0.0);
                    }
                }
            }
        }

        [Test]
        public void DenseOperatorDivide()
        {
            const int dim = 4;
            DiagonalMatrix matrix = new DiagonalMatrix(dim, 4.0);
            matrix = matrix / 2.0;
            for (int row = 0; row < dim; row++)
            {
                for (int col = 0; col < dim; col++)
                {
                    if (row == col)
                    {
                        Assert.AreEqual(2.0, matrix[row, col]);
                    }
                    else
                    {
                        Assert.AreEqual(matrix[row, col], 0.0);
                    }
                }
            }
        }

        [Test]
        public void OperatorAdd()
        {
            const int rows = 4;
            const int cols = 3;
            DenseMatrix matrix1 = new DenseMatrix(rows, cols, 1.0);
            DenseMatrix matrix2 = new DenseMatrix(rows, cols, 2.0);

            DenseMatrix result = matrix1 + matrix2;
            for (int row = 0; row < rows; row++)
            {
                for (int col = 0; col < cols; col++)
                {
                    if (row == col)
                    {
                        Assert.AreEqual(3.0, result[row, col]);
                    }
                }
            }
        }

        [Test]
        public void OperatorAddThrowsExceptions()
        {
            DiagonalMatrix matrix1 = new DiagonalMatrix(2, 2);
            DiagonalMatrix matrix2 = null;
            DiagonalMatrix result;
            Assert.Throws<ArgumentNullException>(() => result = matrix1 + matrix2);
            Assert.Throws<ArgumentNullException>(() => result = matrix2 + matrix1);
            matrix2 = new DiagonalMatrix(3, 3);
            Assert.Throws<NotConformableException>(() => result = matrix1 + matrix2);
        }

        [Test]
        public void OperatorSubtract()
        {
            const int rows = 4;
            const int cols = 3;
            DiagonalMatrix matrix1 = new DiagonalMatrix(rows, cols, 1.0);
            DiagonalMatrix matrix2 = new DiagonalMatrix(rows, cols, 2.0);

            DiagonalMatrix result = matrix1 - matrix2;
            for (int row = 0; row < rows; row++)
            {
                for (int col = 0; col < cols; col++)
                {
                    if (row == col)
                    {
                        Assert.AreEqual(-1.0, result[row, col]);
                    }
                }
            }
        }

        [Test]
        public void OperatorSubtractThrowsExceptions()
        {
            DiagonalMatrix matrix1 = new DiagonalMatrix(2, 2);
            DiagonalMatrix matrix2 = null;
            DiagonalMatrix result;
            Assert.Throws<ArgumentNullException>(() => result = matrix1 - matrix2);
            Assert.Throws<ArgumentNullException>(() => result = matrix2 - matrix1);
            matrix2 = new DiagonalMatrix(3, 3);
            Assert.Throws<NotConformableException>(() => result = matrix1 - matrix2);
        }

        [Test]
        public void OperatorMultiplyScalar()
        {
            DiagonalMatrix m = new DiagonalMatrix(diag3);
            DiagonalMatrix result = .1234 * m;
            for (int i = 0; i < m.Rows; i++)
            {
                for (int j = 0; j < m.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(m[i, j] * .1234, result[i, j], Constants.AcceptableError);
                    }
                }
            }

            m = new DiagonalMatrix(4,3,3.0);
            result = .1234 * m;
            for (int i = 0; i < m.Rows; i++)
            {
                for (int j = 0; j < m.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(m[i, j] * .1234, result[i, j], Constants.AcceptableError);
                    }
                }
            }
        }

        [Test]
        public void OperatorMultiplyMatrix()
        {
            DiagonalMatrix matrix1 = GetMatrix(diag3);
            DiagonalMatrix matrix2 = GetMatrix(diag3);
            Matrix result = matrix1 * matrix2;

            Assert.AreEqual(3, result.Rows);
            Assert.AreEqual(3, result.Columns);
            Assert.AreEqual(1.21, result[0, 0], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[0, 1], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[0, 2], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[1, 0], Constants.AcceptableError);
            Assert.AreEqual(4.84, result[1, 1], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[1, 2], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[2, 0], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[2, 1], Constants.AcceptableError);
            Assert.AreEqual(10.89, result[2, 2], Constants.AcceptableError);
        }
        
        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatotMulitplyMatrixNull()
        {
            DiagonalMatrix matrix1 = null;
            DiagonalMatrix matrix2 = GetMatrix(diag3);
            DiagonalMatrix result = matrix1 * matrix2;
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatotMulitplyMatrixMatrixNull()
        {
            DiagonalMatrix matrix1 = GetMatrix(diag3);
            DiagonalMatrix matrix2 = null;
            DiagonalMatrix result = matrix1 * matrix2;
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void OperatotMulitplyMatrixMatrixNotConformable()
        {
            DiagonalMatrix matrix1 = GetMatrix(diag3);
            DiagonalMatrix matrix2 = new DiagonalMatrix(4,8,1.0);
            DiagonalMatrix result = matrix1 * matrix2;
        }

        [Test]
        public void OperatorMultiplyMatrixVector()
        {
            DiagonalMatrix m = GetMatrix(diag3);
            DenseVector v = new DenseVector(new double[] { 1.1, 2.2, 3.3 });
            Vector result = m * v;
            Assert.AreEqual(v.Count, result.Count);
            Assert.AreEqual(-1.21, result[0], Constants.AcceptableError);
            Assert.AreEqual(-4.84, result[1], Constants.AcceptableError);
            Assert.AreEqual(-10.89, result[2], Constants.AcceptableError);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatorMultiplyMatrixNullVector()
        {
            DiagonalMatrix m = null;
            DenseVector v = new DenseVector(new double[] { 1.1, 2.2, 3.3 });
            Vector result = m * v;
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatorMultiplyMatrixVectorNull()
        {
            DiagonalMatrix m = GetMatrix(diag3);
            DenseVector v = null;
            Vector result = m * v;
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void OperatorMultiplyMatrixVectorNotConformable()
        {
            DiagonalMatrix m = GetMatrix(diag3);
            DenseVector v = new DenseVector(new double[] { 1.1, 2.2 });
            Vector result = m * v;
        }

        [Test]
        public void OperatorMultiplyVector()
        {
            DiagonalMatrix m = GetMatrix(diag3);
            DenseVector v = new DenseVector(new double[] { 1.1, 2.2, 3.3 });
            Vector result = v * m;
            Assert.AreEqual(v.Count, result.Count);
            Assert.AreEqual(-1.21, result[0], Constants.AcceptableError);
            Assert.AreEqual(-4.84, result[1], Constants.AcceptableError);
            Assert.AreEqual(-10.89, result[2], Constants.AcceptableError);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatorMultiplyVectorMatrixNull()
        {
            DiagonalMatrix m = null;
            DenseVector v = new DenseVector(new double[] { 1.1, 2.2, 3.3 });
            Vector result = v * m;
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatorMultiplyVectorNull()
        {
            DiagonalMatrix m = GetMatrix(diag3);
            DenseVector v = null;
            Vector result = v * m;
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void OperatorMultiplyVectorMatrixNotConformable()
        {
            DiagonalMatrix m = GetMatrix(diag3);
            DenseVector v = new DenseVector(new double[] { 1.1, 2.2 });
            Vector result = v * m;
        }

        [Test]
        public void OperatorDivisionThrowsExceptions()
        {
            DenseMatrix matrix = null;
            DenseMatrix result;
            Assert.Throws<ArgumentNullException>(() => result = matrix / 3);
        }

        [TestCase(0)]
        public void AddScalar(double scalar)
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                matrix.Add(scalar);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(GetMatrix(data)[i, j] + scalar, matrix[i, j]);
                    }
                }
            }
        }

        [TestCase(1)]
        [TestCase(2.2)]
        [ExpectedException(typeof(IndexOutOfRangeException))]
        public void AddScalarFail(double scalar)
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                matrix.Add(scalar);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(GetMatrix(data)[i, j] + scalar, matrix[i, j]);
                    }
                }
            }
        }

        [TestCase(0)]
        public void AddScalarResult(double scalar)
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix result = GetMatrix(matrix.Rows, matrix.Columns);
                matrix.Add(scalar, result);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j] + scalar, result[i, j]);
                    }
                }
            }
        }

        [TestCase(1)]
        [TestCase(2.2)]
        [ExpectedException(typeof(IndexOutOfRangeException))]
        public void AddScalarResultFail(double scalar)
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix result = GetMatrix(matrix.Rows, matrix.Columns);
                matrix.Add(scalar, result);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j] + scalar, result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void AddScalarResultNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = null;
            matrix.Add(2.2, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void AddScalarResultNotConformable1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = GetMatrix(matrix.Rows + 1, matrix.Columns);
            matrix.Add(2.2, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void AddScalarResultNotConformable2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = GetMatrix(matrix.Rows, matrix.Columns + 1);
            matrix.Add(2.2, result);
        }

        [Test]
        public void AddMatrix()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix other = matrix.Clone();
                matrix.Add(other);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(other[i, j] + other[i, j], matrix[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void AddMatrixNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = null;
            matrix.Add(other);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void AddMatrixNotConformable1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(matrix.Rows + 1, matrix.Columns);
            matrix.Add(other);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void AddMatrixNotConformable2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(matrix.Rows, matrix.Columns + 1);
            matrix.Add(other);
        }

        [Test]
        [ExpectedException(typeof(ArgumentException))]
        public void AddMatrixNotDiagonal()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix dense = GetNonDiagonalMatrix(3, 3);
            matrix.Add(dense);
        }

        [Test]
        public void AddMatrixResult()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix other = matrix.Clone();
                Matrix result = GetMatrix(matrix.Rows, matrix.Columns);
                matrix.Add(other, result);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j] + other[i, j], result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void AddMatrixResultNull1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = null;
            Matrix result = GetMatrix(matrix.Rows, matrix.Columns);
            matrix.Add(other, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void AddMatrixResultNull2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = matrix.Clone();
            Matrix result = null;
            matrix.Add(other, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void AddMatrixResultNotConformable1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(matrix.Rows + 1, matrix.Columns);
            Matrix result = GetMatrix(matrix.Rows, matrix.Columns);
            matrix.Add(other, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void AddMatrixResultNotConformable2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(matrix.Rows, matrix.Columns + 1);
            Matrix result = GetMatrix(matrix.Rows, matrix.Columns);
            matrix.Add(other, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void AddMatrixResultNotConformable3()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(matrix.Rows, matrix.Columns);
            Matrix result = GetMatrix(matrix.Rows + 1, matrix.Columns);
            matrix.Add(other, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void AddMatrixResultNotConformable4()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(matrix.Rows, matrix.Columns);
            Matrix result = GetMatrix(matrix.Rows, matrix.Columns + 1);
            matrix.Add(other, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentException))]
        public void AddMatrixResultNotDiagonal()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix dense = GetNonDiagonalMatrix(3, 3);
            Matrix result = GetMatrix(3, 3);
            matrix.Add(dense, result);
        }

        [Test]
        public void Plus()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(diag3);
                Matrix result = matrix.Plus();
                Assert.IsFalse(ReferenceEquals(matrix, result));
                Assert.AreEqual(matrix.Rows, result.Rows);
                Assert.AreEqual(matrix.Columns, result.Columns);
                for (int i = 0; i < result.Rows; i++)
                {
                    for (int j = 0; j < result.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j], result[i, j]);
                    }
                }
            }
        }

        [Test]
        public void Clear()
        {
            Matrix matrix = GetMatrix(diag3);
            matrix.Clear();
            for (int i = 0; i < matrix.Rows; i++)
            {
                for (int j = 0; j < matrix.Columns; j++)
                {
                    Assert.AreEqual(0, matrix[i, j]);
                }
            }
        }

        [Test]
        public void Clone()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = matrix.Clone();
            for (int i = 0; i < matrix.Rows; i++)
            {
                for (int j = 0; j < matrix.Columns; j++)
                {
                    Assert.AreEqual(matrix[i, j], other[i, j]);
                }
            }
        }

        [Test]
        public void Append()
        {
            Matrix left = GetMatrix(diag3);
            Matrix right = GetMatrix(singular);
            Matrix result = left.Append(right);
            Assert.AreEqual(left.Columns + right.Columns, result.Columns);
            Assert.AreEqual(left.Rows, right.Rows);

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (j < left.Columns)
                    {
                        Assert.AreEqual(left[i, j], result[i, j]);
                    }
                    else
                    {
                        Assert.AreEqual(right[i, j - left.Columns], result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void AppendNull()
        {
            Matrix left = GetMatrix(diag3);
            Matrix right = null;
            left.Append(right);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void AppendNotConformable()
        {
            Matrix left = GetMatrix(diag7);
            Matrix right = GetMatrix(diag5);
            left.Append(right);
        }

        [Test]
        public void AppendResult()
        {
            Matrix left = GetMatrix(diag3);
            Matrix right = GetMatrix(singular);
            Matrix result = new SparseMatrix(left.Rows, left.Columns + right.Columns);
            left.Append(right, result);
            Assert.AreEqual(left.Columns + right.Columns, result.Columns);
            Assert.AreEqual(left.Rows, right.Rows);

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (j < left.Columns)
                    {
                        Assert.AreEqual(left[i, j], result[i, j]);
                    }
                    else
                    {
                        Assert.AreEqual(right[i, j - left.Columns], result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void AppendResultNull1()
        {
            Matrix left = GetMatrix(diag3);
            Matrix right = null;
            Matrix result = GetMatrix(left.Rows + left.Rows, left.Columns);
            left.Append(right, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void AppendResultNull2()
        {
            Matrix left = GetMatrix(diag3);
            Matrix right = GetMatrix(singular);
            Matrix result = null;
            left.Append(right, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void AppendResultNotConformable1()
        {
            Matrix left = GetMatrix(diag7);
            Matrix right = GetMatrix(diag5);
            Matrix result = GetMatrix(left.Columns + right.Columns, left.Columns);
            left.Append(right, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void AppendResultNotConformable2()
        {
            Matrix left = GetMatrix(singular);
            Matrix right = GetMatrix(diag7);
            Matrix result = GetMatrix(left.Columns + right.Columns + 1, left.Columns);
            left.Append(right, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void AppendResultNotConformable3()
        {
            Matrix left = GetMatrix(diag7);
            Matrix right = GetMatrix(diag7);
            Matrix result = GetMatrix(left.Columns + right.Columns, left.Columns + 1);
            left.Append(right, result);
        }

        [Test]
        public void ICloneable()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = (Matrix)((ICloneable)matrix).Clone();
            for (int i = 0; i < matrix.Rows; i++)
            {
                for (int j = 0; j < matrix.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(diag3[j], matrix[i, j]);
                    }
                    else
                    {
                        Assert.AreEqual(0.0, matrix[i, j]);
                    }
                    Assert.AreEqual(matrix[i, j], other[i, j]);
                }
            }
        }

        [Test]
        public void GetColumn()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                for (int i = 0; i < matrix.Columns; i++)
                {
                    Vector column = matrix.GetColumn(i);
                    Assert.AreEqual(matrix.Rows, column.Count);
                    for (int j = 0; j < matrix.Rows; j++)
                    {
                        Assert.AreEqual(matrix[j, i], column[j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void GetColumnArgumentOutOfRange1()
        {
            Matrix matrix = GetMatrix(diag3);
            matrix.GetColumn(-1);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void GetColumnArgumentOutOfRange2()
        {
            Matrix matrix = GetMatrix(diag3);
            matrix.GetColumn(matrix.Columns);
        }

        [Test]
        public void GetColumnResult()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                for (int i = 0; i < matrix.Columns; i++)
                {
                    Vector column = new DenseVector(matrix.Rows);
                    matrix.GetColumn(i, column);
                    for (int j = 0; j < matrix.Rows; j++)
                    {
                        Assert.AreEqual(matrix[j, i], column[j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void GetColumnResultNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector column = null;
            matrix.GetColumn(0, column);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void GetColumnResultArgumentOutOfRange1()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector column = new DenseVector(matrix.Rows);
            matrix.GetColumn(-1, column);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void GetColumnResultArgumentOutOfRange2()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector column = new DenseVector(matrix.Rows);
            matrix.GetColumn(matrix.Columns, column);
        }

        [TestCase(0, 2)]
        [TestCase(1, 1)]
        [TestCase(2, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(-1, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(4, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, 4, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, -1, ExpectedException = typeof(ArgumentException))]
        public void GetColumnRange(int rowStart, int length)
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                for (int i = 0; i < matrix.Columns; i++)
                {
                    Vector column = matrix.GetColumn(i, rowStart, length);
                    Assert.AreEqual(length, column.Count);
                    for (int j = rowStart, k = 0; j < rowStart + length; j++, k++)
                    {
                        Assert.AreEqual(matrix[j, i], column[k]);
                    }
                }
            }
        }

        [TestCase(0, -1, ExpectedException = typeof(ArgumentException))]
        public void GetColumnRangeResultArgumentExeptions(int rowStart, int length)
        {
            Matrix matrix = GetMatrix(diag3);
            Vector result = new DenseVector(matrix.Rows - 1);
            matrix.GetColumn(0, rowStart, length, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentException))]
        public void GetColumnRangeResultArgumentExeption()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector result = new DenseVector(matrix.Rows);
            matrix.GetColumn(0, 0, 0, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void GetColumnRangeResultNotConformableException()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector result = new DenseVector(matrix.Rows - 1);
            matrix.GetColumn(0, 0, matrix.Columns, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentException))]
        public void GetRowRangeResultArgumentExeption()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector result = new DenseVector(matrix.Columns);
            matrix.GetRow(0, 0, 0, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void GetRowRangeResultNotConformableException()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector result = new DenseVector(matrix.Columns - 1);

            matrix.GetRow(0, 0, matrix.Rows, result);
        }

        [Test]
        public void ColumnEnumerator()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                foreach (KeyValuePair<int, Vector> column in matrix.GetColumnEnumerator())
                {
                    int i = 0;
                    foreach (double value in column.Value)
                    {
                        Assert.AreEqual(matrix[i++, column.Key], value);
                    }
                }
            }
        }

        [TestCase(0, 2)]
        [TestCase(1, 1)]
        [TestCase(-1, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(4, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, 4, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, -1, ExpectedException = typeof(ArgumentException))]
        public void ColumnEnumeratorRange(int colStart, int count)
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                foreach (KeyValuePair<int, Vector> column in matrix.GetColumnEnumerator(colStart, count))
                {
                    int i = 0;
                    foreach (double value in column.Value)
                    {
                        Assert.AreEqual(matrix[i++, column.Key], value);
                    }
                }
            }
        }

        [TestCase(3, 4)]
        [TestCase(4, 4)]
        [TestCase(4, 3)]
        public void RowsColumns(int rows, int columns)
        {
            Matrix matrix = GetMatrix(rows, columns);
            Assert.AreEqual(rows, matrix.Rows);
            Assert.AreEqual(columns, matrix.Columns);
        }

        [Test]
        public void ConditionNumber()
        {
            Matrix matrix = GetMatrix(diag3);
            double conditionNumber = matrix.ConditionNumber();
            Assert.AreEqual(3.3 / 1.1, conditionNumber, Constants.AcceptableError);

            matrix = GetMatrix(diag7);
            conditionNumber = matrix.ConditionNumber();
            Assert.AreEqual(7.7 / 1.1, conditionNumber, Constants.AcceptableError);

            matrix = GetMatrix(singular);
            conditionNumber = matrix.ConditionNumber();
            Assert.AreEqual(System.Double.PositiveInfinity, conditionNumber, Constants.AcceptableError);
        }


        [Test]
        public void CopyTo()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix target = GetMatrix(matrix.Rows, matrix.Columns);
                matrix.CopyTo(target);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j], target[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void CopyToNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix target = null;
            matrix.CopyTo(target);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void CopyToNotConformableException()
        {
            Matrix matrix = GetMatrix(diag7);
            Matrix target = GetMatrix(matrix.Columns, matrix.Rows+1);
            matrix.CopyTo(target);
        }

        [Test]
        public void Determinant()
        {
            Matrix matrix = GetMatrix(diag3);
            double det = matrix.Determinant();
            Assert.AreEqual(-1.1*2.2*3.3, det, Constants.AcceptableError);
        }

        [Test]
        [ExpectedException(typeof(MatrixNotSquareException))]
        public void DeterminantNotSquareMatrix()
        {
            Matrix matrix = new DiagonalMatrix(3, 4, 5.0);
            double det = matrix.Determinant();
        }

        [Test]
        public void Trace()
        {
            Matrix matrix = GetMatrix(diag3);
            double trace = matrix.Trace();
            Assert.AreEqual(-1.1-2.2-3.3, trace, Constants.AcceptableError);
        }

        [Test]
        [ExpectedException(typeof(MatrixNotSquareException))]
        public void TracNotSquareMatrix()
        {
            Matrix matrix = new DiagonalMatrix(3, 4, 5.0);
            double trace = matrix.Trace();
        }

        [Test]
        public void Diagonal()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Vector vector = matrix.GetDiagonal();
                int size = System.Math.Min(matrix.Rows, matrix.Columns);
                Assert.AreEqual(size, vector.Count);
                for (int i = 0; i < size; i++)
                {
                    Assert.AreEqual(matrix[i, i], vector[i]);
                }
            }
        }

        [Test]
        public void DiagonalStack()
        {
            Matrix top = GetMatrix(diag5);
            Matrix bottom = GetMatrix(diag7);
            Matrix result = top.DiagonalStack(bottom);
            Assert.AreEqual(top.Rows + bottom.Rows, result.Rows);
            Assert.AreEqual(top.Columns + bottom.Columns, result.Columns);

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i < top.Rows && j < top.Columns)
                    {
                        Assert.AreEqual(top[i, j], result[i, j]);
                    }
                    else if (i >= top.Rows && j >= top.Columns)
                    {
                        Assert.AreEqual(bottom[i - top.Rows, j - top.Columns], result[i, j]);
                    }
                    else
                    {
                        Assert.AreEqual(0, result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void DiagonalStackNull()
        {
            Matrix top = GetMatrix(diag3);
            Matrix bottom = null;
            top.DiagonalStack(bottom);
        }

        [Test]
        public void DiagonalStackResult()
        {
            Matrix top = GetMatrix(diag5);
            Matrix bottom = GetMatrix(diag7);
            Matrix result = GetMatrix(top.Rows + bottom.Rows, top.Columns + bottom.Columns);
            top.DiagonalStack(bottom, result);

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i < top.Rows && j < top.Columns)
                    {
                        Assert.AreEqual(top[i, j], result[i, j]);
                    }
                    else if (i >= top.Rows && j >= top.Columns)
                    {
                        Assert.AreEqual(bottom[i - top.Rows, j - top.Columns], result[i, j]);
                    }
                    else
                    {
                        Assert.AreEqual(0, result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void DiagonalStackResultNull1()
        {
            Matrix top = GetMatrix(diag3);
            Matrix bottom = null;
            Matrix result = GetMatrix(top.Rows + top.Rows, top.Columns + top.Columns);
            top.DiagonalStack(bottom, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void DiagonalStackResultNull2()
        {
            Matrix top = GetMatrix(diag3);
            Matrix bottom = GetMatrix(diag3);
            Matrix result = null;
            top.DiagonalStack(bottom, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void DiagonalStackResultNotConformableException1()
        {
            Matrix top = GetMatrix(diag3);
            Matrix bottom = GetMatrix(diag3);
            
            Matrix result = GetMatrix(top.Rows + bottom.Rows + 1, top.Columns + bottom.Columns);
            top.DiagonalStack(bottom, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void DiagonalStackResultNotConformableException2()
        {
            Matrix top = GetMatrix(diag3);
            Matrix bottom = GetMatrix(diag3);
            
            Matrix result = GetMatrix(top.Rows + bottom.Rows, top.Columns + bottom.Columns + 1);
            top.DiagonalStack(bottom, result);
        }

        [Test]
        public void Divide()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                matrix.Divide(2.3);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        if (i == j)
                        {
                            Assert.AreEqual(data[i] / 2.3, matrix[i, j], Constants.AcceptableError);
                        }
                        else
                        {
                            Assert.AreEqual(0.0, matrix[i, j], Constants.AcceptableError);
                        }
                    }
                }
            }
        }

        [Test]
        public void DivideResult()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix result = GetMatrix(matrix.Rows, matrix.Columns);
                matrix.Divide(2.3, result);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        if (i == j)
                        {
                            Assert.AreEqual(data[i] / 2.3, result[i, j], Constants.AcceptableError);
                        }
                        else
                        {
                            Assert.AreEqual(0.0, result[i, j], Constants.AcceptableError);
                        }
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void DivideResultNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = null;
            matrix.Divide(2.3, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void DivideResultNotConformable1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = GetMatrix(matrix.Rows + 1, matrix.Columns);
            matrix.Divide(2.3, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void DivideResultNotConformable2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = GetMatrix(matrix.Rows, matrix.Columns + 1);
            matrix.Divide(2.3, result);
        }

        [Test]
        public void FrobeniusNorm()
        {
            Matrix matrix = GetMatrix(diag3);
            Assert.AreEqual(System.Math.Sqrt((1.1*1.1) + (2.2*2.2) + (3.3*3.3)), matrix.FrobeniusNorm(), Constants.AcceptableError);

            matrix = GetMatrix(diag5);
            Assert.AreEqual(System.Math.Sqrt((1.1 * 1.1) + (2.2 * 2.2) + (3.3 * 3.3) + (4.4 * 4.4) + (5.5 * 5.5)),
                matrix.FrobeniusNorm(), Constants.AcceptableError);
        }

        [Test]
        public void InfinityNorm()
        {
            Matrix matrix = GetMatrix(diag3);
            Assert.AreEqual(3.3, matrix.InfinityNorm());

            matrix = GetMatrix(diag5);
            Assert.AreEqual(5.5, matrix.InfinityNorm());

            matrix = GetMatrix(diag7);
            Assert.AreEqual(7.7, matrix.InfinityNorm());
        }

        [Test]
        public void L1Norm()
        {
            Matrix matrix = GetMatrix(diag3);
            Assert.AreEqual(3.3, matrix.L1Norm());

            matrix = GetMatrix(diag5);
            Assert.AreEqual(5.5, matrix.L1Norm());

            matrix = GetMatrix(diag7);
            Assert.AreEqual(7.7, matrix.L1Norm());
        }

        [Test]
        public void L2Norm()
        {
            Matrix matrix = GetMatrix(diag3);
            Assert.AreEqual(3.3, matrix.L2Norm(), Constants.AcceptableError);

            matrix = GetMatrix(diag5);
            Assert.AreEqual(5.5, matrix.L2Norm(), Constants.AcceptableError);

            matrix = GetMatrix(diag7);
            Assert.AreEqual(7.7, matrix.L2Norm(), Constants.AcceptableError);
        }

        [Test]
        public void Transpose()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix result = matrix.Transpose();

                Assert.AreEqual(matrix.Rows, result.Columns);
                Assert.AreEqual(matrix.Columns, result.Rows);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j], result[j, i]);
                    }
                }
            }
        }

        [Test]
        public void Inverse()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix inverse = matrix.Inverse();

            Assert.AreEqual(matrix.Rows, inverse.Rows);
            Assert.AreEqual(matrix.Columns, inverse.Columns);
            Matrix identity = new DenseMatrix(matrix.Rows, matrix.Columns);
            matrix.Multiply(inverse, identity);

            for (int i = 0; i < identity.Rows; i++)
            {
                for (int j = 0; j < identity.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(1.0, identity[i, j], Constants.AcceptableError);
                    }
                    else
                    {
                        Assert.AreEqual(0.0, identity[i, j], Constants.AcceptableError);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void InverseSingular()
        {
            Matrix singularMatrix = GetMatrix(singular);
            singularMatrix.Inverse();
        }

        [Test]
        [ExpectedException(typeof(MatrixNotSquareException))]
        public void InverseNotSquareMatrix()
        {
            Matrix matrix = GetNonDiagonalMatrix(4, 2);
            matrix.Inverse();
        }

        [Test]
        public void LeftMultiplyVector()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector vector = GetVector(new double[] { 1.1, 2.2, 3.3 });
            Vector result = matrix.LeftMultiply(vector);
            Assert.AreEqual(matrix.Columns, result.Count);
            Assert.AreEqual(-1.1*1.1, result[0], Constants.AcceptableError);
            Assert.AreEqual(-2.2*2.2, result[1], Constants.AcceptableError);
            Assert.AreEqual(-3.3*3.3, result[2], Constants.AcceptableError);

            matrix = GetMatrix(singular);
            vector = GetVector(new double[] { 1.1, 2.2, 3.3 });
            result = matrix.LeftMultiply(vector);
            Assert.AreEqual(matrix.Columns, result.Count);
            Assert.AreEqual(-3.3*1.1, result[0], Constants.AcceptableError);
            Assert.AreEqual(-4.4*2.2, result[1], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[2], Constants.AcceptableError);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void LeftMultiplyVectorNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector leftSide = null;
            matrix.LeftMultiply(leftSide);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void LeftMultiplyVectorNotConformable()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector leftSide = GetVector(matrix.Rows + 1);
            matrix.LeftMultiply(leftSide);
        }

        [Test]
        public void LeftMultiplyVectorResult()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector vector = GetVector(new double[] { 1.1, 2.2, 3.3 });
            Vector result = GetVector(matrix.Columns);
            matrix.LeftMultiply(vector, result);
            Assert.AreEqual(-1.1 * 1.1, result[0], Constants.AcceptableError);
            Assert.AreEqual(-2.2 * 2.2, result[1], Constants.AcceptableError);
            Assert.AreEqual(-3.3 * 3.3, result[2], Constants.AcceptableError);

            matrix = GetMatrix(singular);
            vector = GetVector(new double[] { 1.1, 2.2, 3.3 });
            matrix.LeftMultiply(vector, result);
            Assert.AreEqual(matrix.Columns, result.Count);
            Assert.AreEqual(-3.3 * 1.1, result[0], Constants.AcceptableError);
            Assert.AreEqual(-4.4 * 2.2, result[1], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[2], Constants.AcceptableError);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void LeftMultiplyVectorResultNull1()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector leftSide = null;
            Vector result = GetVector(matrix.Columns);
            matrix.LeftMultiply(leftSide, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void LeftMultiplyVectorResultNull2()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector leftSide = GetVector(matrix.Rows);
            Vector result = null;
            matrix.LeftMultiply(leftSide, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void LeftMultiplyVectorResultNotConformable1()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector leftSide = GetVector(matrix.Rows + 1);
            
            Vector result = GetVector(matrix.Columns);
            matrix.LeftMultiply(leftSide, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void LeftMultiplyVectorResultNotConformable2()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector leftSide = GetVector(matrix.Rows);
            ;
            Vector result = GetVector(matrix.Columns + 1);
            matrix.LeftMultiply(leftSide, result);
        }

        [Test]
        public void LowerTriangle()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix lower = matrix.GetLowerTriangle();
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        if (i >= j)
                        {
                            Assert.AreEqual(matrix[i, j],
                                            lower[i, j]);
                        }
                        else
                        {
                            Assert.AreEqual(0, lower[i, j]);
                        }
                    }
                }
            }
        }


        [Test]
        public void LowerTriangleResult()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix lower = GetMatrix(matrix.Rows, matrix.Columns);
                matrix.GetLowerTriangle(lower);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        if (i >= j)
                        {
                            Assert.AreEqual(matrix[i, j], lower[i, j]);
                        }
                        else
                        {
                            Assert.AreEqual(0, lower[i, j]);
                        }
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void LowerTriangleResultNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix lower = null;
            matrix.GetLowerTriangle(lower);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void LowerTriangleResultNotConformableException1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix lower = GetMatrix(matrix.Rows + 1, matrix.Columns);
            matrix.GetLowerTriangle(lower);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void LowerTriangleResultNotConformableException2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix lower = GetMatrix(matrix.Rows, matrix.Columns + 1);
            matrix.GetLowerTriangle(lower);
        }

        [Test]
        public void StrictlyLowerTriangle()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix lower = matrix.GetStrictlyLowerTriangle();
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        if (i > j)
                        {
                            Assert.AreEqual(matrix[i, j], lower[i, j]);
                        }
                        else
                        {
                            Assert.AreEqual(0, lower[i, j]);
                        }
                    }
                }
            }
        }

        [Test]
        public void StrictlyLowerTriangleResult()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix lower = GetMatrix(matrix.Rows, matrix.Columns);
                matrix.GetStrictlyLowerTriangle(lower);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        if (i > j)
                        {
                            Assert.AreEqual(matrix[i, j], lower[i, j]);
                        }
                        else
                        {
                            Assert.AreEqual(0, lower[i, j]);
                        }
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void StrictlyLowerTriangleResultNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix lower = null;
            matrix.GetStrictlyLowerTriangle(lower);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void StrictlyLowerTriangleResultNotConformableException1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix lower = GetMatrix(matrix.Rows + 1, matrix.Columns);
            matrix.GetStrictlyLowerTriangle(lower);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void StrictlyLowerTriangleResultNotConformableException2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix lower = GetMatrix(matrix.Rows, matrix.Columns + 1);
            matrix.GetStrictlyLowerTriangle(lower);
        }

        [Test]
        public void UpperTriangle()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix upper = matrix.GetUpperTriangle();
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        if (i <= j)
                        {
                            Assert.AreEqual(matrix[i, j], upper[i, j]);
                        }
                        else
                        {
                            Assert.AreEqual(0, upper[i, j]);
                        }
                    }
                }
            }
        }

        [Test]
        public void UpperTriangleResult()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix upper = GetMatrix(matrix.Rows, matrix.Columns);
                matrix.GetUpperTriangle(upper);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        if (i <= j)
                        {
                            Assert.AreEqual(matrix[i, j], upper[i, j]);
                        }
                        else
                        {
                            Assert.AreEqual(0, upper[i, j]);
                        }
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void UpperTriangleResultNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix upper = null;
            matrix.GetUpperTriangle(upper);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void UpperTriangleResultNotConformableException1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix upper = GetMatrix(matrix.Rows + 1, matrix.Columns);
            matrix.GetUpperTriangle(upper);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void UpperTriangleResultNotConformableException2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix upper = GetMatrix(matrix.Rows, matrix.Columns + 1);
            matrix.GetUpperTriangle(upper);
        }

        [Test]
        public void StrictlyUpperTriangle()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix upper = matrix.GetStrictlyUpperTriangle();
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        if (i < j)
                        {
                            Assert.AreEqual(matrix[i, j], upper[i, j]);
                        }
                        else
                        {
                            Assert.AreEqual(0, upper[i, j]);
                        }
                    }
                }
            }
        }

        [Test]
        public void StrictlyUpperTriangleResult()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix upper = GetMatrix(matrix.Rows, matrix.Columns);
                matrix.GetStrictlyUpperTriangle(upper);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        if (i < j)
                        {
                            Assert.AreEqual(matrix[i, j], upper[i, j]);
                        }
                        else
                        {
                            Assert.AreEqual(0, upper[i, j]);
                        }
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void StrictlyUpperTriangleResultNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix upper = null;
            matrix.GetStrictlyUpperTriangle(upper);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void StrictlyUpperTriangleResultNotConformableException1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix upper = GetMatrix(matrix.Rows + 1, matrix.Columns);
            matrix.GetStrictlyUpperTriangle(upper);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void StrictlyUpperTriangleResultNotConformableException2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix upper = GetMatrix(matrix.Rows, matrix.Columns + 1);
            matrix.GetStrictlyUpperTriangle(upper);
        }

        [TestCase(0)]
        [TestCase(1)]
        [TestCase(2.2)]
        public void MultiplyScalar(double scalar)
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix matrixs = GetMatrix(data);
                matrixs.Multiply(scalar);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j] * scalar, matrixs[i, j], Constants.AcceptableError);
                    }
                }
            }
        }

        [TestCase(0)]
        [TestCase(1)]
        [TestCase(2.2)]
        public void MultiplyScalarResult(double scalar)
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix result = GetMatrix(matrix.Rows, matrix.Columns);
                matrix.Multiply(scalar, result);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j] * scalar, result[i, j], Constants.AcceptableError);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void MultiplyScalarResultNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = null;
            matrix.Multiply(2.3, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void MultiplyScalarResultNotConformable1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = GetMatrix(matrix.Rows + 1, matrix.Columns);
            matrix.Multiply(2.3, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void MultiplyScalarResultNotConformable2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = GetMatrix(matrix.Rows, matrix.Columns + 1);
            matrix.Multiply(2.3, result);
        }

        [Test]
        public void MultiplyMatrix()
        {
            Matrix matrix1 = GetMatrix(diag3);
            Matrix matrix2 = matrix1.Clone();
            Matrix result = matrix1.Multiply(matrix2);

            Assert.AreEqual(3, result.Rows);
            Assert.AreEqual(3, result.Columns);
            Assert.AreEqual(matrix1[0, 0] * matrix2[0, 0], result[0, 0], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[0, 1], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[0, 2], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[1, 0], Constants.AcceptableError);
            Assert.AreEqual(matrix1[1, 1] * matrix2[1, 1], result[1, 1], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[1, 2], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[2, 0], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[2, 1], Constants.AcceptableError);
            Assert.AreEqual(matrix1[2, 2] * matrix2[2, 2], result[2, 2], Constants.AcceptableError);

            matrix2 = new DenseMatrix(new double[,] { {1.0,2.0,3.0,4.0}, {5.0, 6.0, 7.0, 8.0}, {9.0, 10.0, 11.0, 12.0}});
            result = matrix1.Multiply(matrix2);
            Assert.AreEqual(3, result.Rows);
            Assert.AreEqual(4, result.Columns);
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    Assert.AreEqual(result[i,j], matrix2[i,j] * matrix1[i,i], Constants.AcceptableError);
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void MultiplyMatrixNull()
        {
            Matrix matrix1 = GetMatrix(diag3);
            Matrix matrix2 = null;
            matrix1.Multiply(matrix2);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void MultiplyMatrixNotConformable()
        {
            Matrix matrix1 = GetMatrix(diag3);
            Matrix matrix2 = GetMatrix(matrix1.Rows + 1, matrix1.Columns);
            matrix1.Multiply(matrix2);
        }

        [Test]
        public void MultiplyMatrixResult()
        {
            Matrix matrix1 = GetMatrix(diag3);
            Matrix matrix2 = GetMatrix(diag3);
            Matrix result = GetMatrix(matrix1.Rows, matrix2.Columns);
            matrix1.Multiply(matrix2, result);
            Assert.AreEqual(3, result.Rows);
            Assert.AreEqual(3, result.Columns);
            Assert.AreEqual(matrix1[0, 0] * matrix2[0, 0], result[0, 0], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[0, 1], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[0, 2], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[1, 0], Constants.AcceptableError);
            Assert.AreEqual(matrix1[1, 1] * matrix2[1, 1], result[1, 1], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[1, 2], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[2, 0], Constants.AcceptableError);
            Assert.AreEqual(0.0, result[2, 1], Constants.AcceptableError);
            Assert.AreEqual(matrix1[2, 2] * matrix2[2, 2], result[2, 2], Constants.AcceptableError);

            result = new DenseMatrix(3, 4);
            matrix2 = new DenseMatrix(new double[,] { { 1.0, 2.0, 3.0, 4.0 }, { 5.0, 6.0, 7.0, 8.0 }, { 9.0, 10.0, 11.0, 12.0 } });
            matrix1.Multiply(matrix2, result);
            Assert.AreEqual(3, result.Rows);
            Assert.AreEqual(4, result.Columns);
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    Assert.AreEqual(result[i, j], matrix2[i, j] * matrix1[i, i], Constants.AcceptableError);
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void MultiplyMatrixResultNull1()
        {
            Matrix matrix1 = GetMatrix(diag3);
            Matrix matrix2 = null;
            Matrix result = GetMatrix(matrix1.Rows, matrix1.Columns);
            matrix1.Multiply(matrix2, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void MultiplyMatrixResultNull2()
        {
            Matrix matrix1 = GetMatrix(diag3);
            Matrix matrix2 = GetMatrix(diag3);
            Matrix result = null;
            matrix1.Multiply(matrix2, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void MultiplyMatrixResultNotConformable1()
        {
            Matrix matrix1 = GetMatrix(diag3);
            Matrix matrix2 = GetMatrix(matrix1.Rows + 1, matrix1.Columns);
            Matrix result = GetMatrix(matrix1.Rows, matrix2.Columns);
            matrix1.Multiply(matrix2, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void MultiplyMatrixResultNotConformable2()
        {
            Matrix matrix1 = GetMatrix(diag3);
            Matrix matrix2 = GetMatrix(matrix1.Rows, matrix1.Columns);
            Matrix result = GetMatrix(matrix1.Rows + 1, matrix2.Columns);
            matrix1.Multiply(matrix2, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void MultiplyMatrixResultNotConformable3()
        {
            Matrix matrix1 = GetMatrix(diag3);
            Matrix matrix2 = GetMatrix(matrix1.Rows, matrix1.Columns);
            Matrix result = GetMatrix(matrix1.Rows, matrix2.Columns + 1);
            matrix1.Multiply(matrix2, result);
        }

        [Test]
        public void MultiplyVector()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector vector = GetVector(new double[] { 1.1, 2.2, 3.3 });
            Vector result = matrix.Multiply(vector);
            Assert.AreEqual(matrix.Rows, result.Count);
            Assert.AreEqual(-1.1*1.1, result[0], Constants.AcceptableError);
            Assert.AreEqual(-2.2*2.2, result[1], Constants.AcceptableError);
            Assert.AreEqual(-3.3*3.3, result[2], Constants.AcceptableError);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void MultiplyVectorNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector rightSide = null;
            matrix.Multiply(rightSide);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void MultiplyVectorNotConformable()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector rightSide = GetVector(matrix.Columns + 1);
            matrix.Multiply(rightSide);
        }

        [Test]
        public void MultiplyVectorResult()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector vector = GetVector(new double[] { 1.1, 2.2, 3.3 });
            Vector result = GetVector(matrix.Rows);
            matrix.Multiply(vector, result);
            Assert.AreEqual(matrix.Rows, result.Count);
            Assert.AreEqual(-1.1 * 1.1, result[0], Constants.AcceptableError);
            Assert.AreEqual(-2.2 * 2.2, result[1], Constants.AcceptableError);
            Assert.AreEqual(-3.3 * 3.3, result[2], Constants.AcceptableError);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void MultiplyVectorResultNull1()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector rightSide = null;
            Vector result = GetVector(matrix.Rows);
            matrix.Multiply(rightSide, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void MultiplyVectorResultNull2()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector rightSide = GetVector(matrix.Columns);
            Vector result = null;
            matrix.Multiply(rightSide, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void MultiplyVectorResultNotConformable1()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector rightSide = GetVector(matrix.Columns + 1);
            Vector result = GetVector(matrix.Rows);
            matrix.Multiply(rightSide, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void MultiplyVectorResultNotConformable2()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector rightSide = GetVector(matrix.Columns);
            Vector result = GetVector(matrix.Rows + 1);
            matrix.Multiply(rightSide, result);
        }

        [Test]
        public void Negate()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                matrix.Negate();
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        if (i == j)
                        {
                            Assert.AreEqual(-data[i], matrix[i, j]);
                        }
                        else
                        {
                            Assert.AreEqual(0.0, matrix[i, j]);
                        }
                    }
                }
            }
        }

        [Test]
        public void NegateResult()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix result = GetMatrix(matrix.Rows, matrix.Columns);
                matrix.Negate(result);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(-matrix[i, j], result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void NegateResultNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = null;
            matrix.Negate(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void NegateResultNotConformable1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = GetMatrix(matrix.Rows + 1, matrix.Columns);
            matrix.Negate(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void NegateResultNotConformable2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = GetMatrix(matrix.Rows, matrix.Columns + 1);
            matrix.Negate(result);
        }

        [Test]
        public void Row()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    Vector row = matrix.GetRow(i);

                    Assert.AreEqual(matrix.Columns, row.Count);
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j], row[j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void RowArgumentOutOfRange1()
        {
            Matrix matrix = GetMatrix(diag3);
            matrix.GetRow(-1);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void RowArgumentOutOfRange2()
        {
            Matrix matrix = GetMatrix(diag3);
            matrix.GetRow(matrix.Columns);
        }

        [Test]
        public void RowResult()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    Vector row = GetVector(matrix.Columns);
                    matrix.GetRow(i, row);
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j], row[j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void RowResultNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector row = null;
            matrix.GetRow(0, row);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void RowResultArgumentOutOfRange1()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector row = GetVector(matrix.Columns);
            matrix.GetRow(-1, row);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void RowResultArgumentOutOfRange2()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector row = GetVector(matrix.Columns);
            matrix.GetRow(matrix.Columns, row);
        }


        [TestCase(0, 2)]
        [TestCase(1, 1)]
        [TestCase(-1, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(4, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, 4, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, -1, ExpectedException = typeof(ArgumentException))]
        public void RowRange(int start, int length)
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    Vector row = matrix.GetRow(i, start, length);
                    Assert.AreEqual(length, row.Count);
                    for (int j = start, k = 0; j < length; j++, k++)
                    {
                        Assert.AreEqual(matrix[i, j], row[k]);
                    }
                }
            }
        }


        [Test]
        public void RowEnumerator()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                foreach (KeyValuePair<int, Vector> row in matrix.GetRowEnumerator())
                {
                    int j = 0;
                    foreach (double value in row.Value)
                    {
                        Assert.AreEqual(matrix[row.Key, j++], value);
                    }
                }
            }
        }

        [TestCase(0, 2)]
        [TestCase(1, 1)]
        [TestCase(-1, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(4, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, 4, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, -1, ExpectedException = typeof(ArgumentException))]
        public void RowEnumeratorRange(int start, int length)
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);

                foreach (KeyValuePair<int, Vector> row in matrix.GetRowEnumerator(start, length))
                {
                    int j = 0;
                    foreach (double value in row.Value)
                    {
                        Assert.AreEqual(matrix[row.Key, j++], value);
                    }
                }
            }
        }

        [Test]
        public void SetColumnVector()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 0.0, 0.0, 0.0 };
            for (int i = 0; i < matrix.Columns; i++)
            {
                array[i] = 1.0;
                Vector vector = GetVector(array);
                matrix.SetColumn(i, vector);
                for (int j = 0; j < matrix.Rows; j++)
                {
                    Assert.AreEqual(array[j], matrix[j, i]);
                }
                array[i] = 0.0;
            }
        }

        [Test]
        [ExpectedException(typeof(IndexOutOfRangeException))]
        public void SetColumnVectorNotDiagonal()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1.0, 1.0, 1.0 };
            Vector vector = GetVector(array);
            matrix.SetColumn(0, vector);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SetColumnVectorNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector vector = null;
            matrix.SetColumn(0, vector);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void SetColumnVectorArgumentOutOfRange1()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1, 2, 3 };
            Vector vector = GetVector(array);
            matrix.SetColumn(-1, vector);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void SetColumnVectorArgumentOutOfRange2()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1, 2, 3 };
            Vector vector = GetVector(array);
            matrix.SetColumn(matrix.Columns + 1, vector);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SetColumnVectorNotConformableException()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1, 2, 3, 4 };
            Vector vector = GetVector(array);
            matrix.SetColumn(0, vector);
        }

        [Test]
        public void SetColumnArray()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 0.0, 0.0, 0.0 };
            for (int i = 0; i < matrix.Columns; i++)
            {
                array[i] = 1.0;
                matrix.SetColumn(i, array);
                for (int j = 0; j < matrix.Rows; j++)
                {
                    Assert.AreEqual(array[j], matrix[j, i]);
                }
                array[i] = 0.0;
            }
        }

        [Test]
        [ExpectedException(typeof(IndexOutOfRangeException))]
        public void SetColumnArrayNotDiagonal()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1.0, 1.0, 1.0 };
            matrix.SetColumn(0, array);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SetColumnArrayNull()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = null;
            matrix.SetColumn(0, array);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void SetColumnArrayArgumentOutOfRange1()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1, 2, 3 };
            matrix.SetColumn(-1, array);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void SetColumnArrayArgumentOutOfRange2()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1, 2, 3 };
            matrix.SetColumn(matrix.Columns + 1, array);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SetColumnArrayNotConformableException()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1, 2, 3, 4 };
            matrix.SetColumn(0, array);
        }

        [Test]
        public void SetRowVector()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 0.0, 0.0, 0.0 };
            for (int i = 0; i < matrix.Rows; i++)
            {
                array[i] = 1.0;
                Vector vector = GetVector(array);
                matrix.SetRow(i, vector);
                for (int j = 0; j < matrix.Columns; j++)
                {
                    Assert.AreEqual(array[j], matrix[i, j]);
                }
                array[i] = 0.0;
            }
        }

        [Test]
        [ExpectedException(typeof(IndexOutOfRangeException))]
        public void SetRowVectorNotDiagonal()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1.0, 1.0, 1.0 };
            Vector vector = GetVector(array);
            matrix.SetRow(0, vector);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SetRowVectorNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector vector = null;
            matrix.SetRow(0, vector);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void SetRowVectorArgumentOutOfRange1()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1, 2, 3 };
            Vector vector = GetVector(array);
            matrix.SetRow(-1, vector);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void SetRowVectorArgumentOutOfRange2()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1, 2, 3 };
            Vector vector = GetVector(array);
            matrix.SetRow(matrix.Rows + 1, vector);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SetRowVectorNotConformableException()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1, 2, 3, 4 };
            Vector vector = GetVector(array);
            matrix.SetRow(0, vector);
        }

        [Test]
        public void SetRowArray()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 0.0, 0.0, 0.0 };
            for (int i = 0; i < matrix.Rows; i++)
            {
                array[i] = 1.0;
                matrix.SetRow(i, array);
                for (int j = 0; j < matrix.Columns; j++)
                {
                    Assert.AreEqual(array[j], matrix[i, j]);
                }
                array[i] = 0.0;
            }
        }

        [Test]
        [ExpectedException(typeof(IndexOutOfRangeException))]
        public void SetRowArrayNotDiagonal()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1.0, 1.0, 1.0 };
            matrix.SetRow(0, array);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SetRowArrayNull()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = null;
            matrix.SetRow(0, array);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void SetRowArrayArgumentOutOfRange1()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1, 2, 3 };
            matrix.SetRow(-1, array);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void SetRowArrayArgumentOutOfRange2()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1, 2, 3 };
            matrix.SetRow(matrix.Rows + 1, array);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SetRowArrayNotConformableException()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1, 2, 3, 4 };
            matrix.SetRow(0, array);
        }

        [Test]
        public void SetDiagonalVector()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1, 2, 3 };
            Vector vector = GetVector(array);
            matrix.SetDiagonal(vector);
            for (int i = 0; i < array.Length; i++)
            {
                Assert.AreEqual(array[i], matrix[i, i]);
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SetDiagonalVectorNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector vector = null;
            matrix.SetDiagonal(vector);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SetDiagonalVectorNotConformable()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector vector = GetVector(matrix.Rows + 1);
            matrix.SetDiagonal(vector);
        }

        [Test]
        public void SetDiagonalArray()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[] { 1, 2, 3 };
            matrix.SetDiagonal(array);
            for (int i = 0; i < array.Length; i++)
            {
                Assert.AreEqual(array[i], matrix[i, i]);
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SetDiagonalArrayNull()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = null;
            matrix.SetDiagonal(array);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SetDiagonalArrayNotConformable()
        {
            Matrix matrix = GetMatrix(diag3);
            double[] array = new double[matrix.Rows + 1];
            matrix.SetDiagonal(array);
        }

        [Test]
        public void Stack()
        {
            Matrix top = GetMatrix(diag3);
            Matrix bottom = GetMatrix(singular);
            Matrix result = top.Stack(bottom);
            Assert.AreEqual(top.Rows + bottom.Rows, result.Rows);
            Assert.AreEqual(top.Columns, result.Columns);

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i < top.Rows)
                    {
                        Assert.AreEqual(top[i, j], result[i, j]);
                    }
                    else
                    {
                        Assert.AreEqual(bottom[i - top.Rows, j], result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void StackNull()
        {
            Matrix top = GetMatrix(diag3);
            Matrix bottom = null;
            top.Stack(bottom);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void StackNotConformable()
        {
            Matrix top = GetMatrix(diag7);
            Matrix bottom = GetMatrix(diag5);
            top.Stack(bottom);
        }

        [Test]
        public void StackResult()
        {
            Matrix top = GetMatrix(diag3);
            Matrix bottom = GetMatrix(singular);
            Matrix result = GetNonDiagonalMatrix(top.Rows + bottom.Rows, top.Columns);
            top.Stack(bottom, result);
            Assert.AreEqual(top.Rows + bottom.Rows, result.Rows);
            Assert.AreEqual(top.Columns, result.Columns);

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i < top.Rows)
                    {
                        Assert.AreEqual(top[i, j], result[i, j]);
                    }
                    else
                    {
                        Assert.AreEqual(bottom[i - top.Rows, j], result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void StackResultNull1()
        {
            Matrix top = GetMatrix(diag3);
            Matrix bottom = null;
            Matrix result = GetMatrix(top.Rows + top.Rows, top.Columns);
            top.Stack(bottom, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void StackResultNull2()
        {
            Matrix top = GetMatrix(diag3);
            Matrix bottom = GetMatrix(diag7);
            Matrix result = null;
            top.Stack(bottom, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void StackResultNotConformable1()
        {
            Matrix top = GetMatrix(diag7);
            Matrix bottom = GetMatrix(diag5);
            Matrix result = GetMatrix(top.Rows + bottom.Rows, top.Columns);
            top.Stack(bottom, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void StackResultNotConformable2()
        {
            Matrix top = GetMatrix(diag7);
            Matrix bottom = GetMatrix(diag7);
            Matrix result = GetMatrix(top.Rows + bottom.Rows + 1, top.Columns);
            top.Stack(bottom, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void StackResultNotConformable3()
        {
            Matrix top = GetMatrix(diag7);
            Matrix bottom = GetMatrix(diag7);
            Matrix result = GetMatrix(top.Rows + bottom.Rows, top.Columns + 1);
            top.Stack(bottom, result);
        }

        [TestCase(0, 2, 0, 2)]
        [TestCase(1, 1, 1, 1)]
        [TestCase(0, 4, 0, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, 2, 0, 4, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(4, 2, 0, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, 2, 4, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(-1, 2, 0, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, 2, -1, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, -1, 0, 2, ExpectedException = typeof(ArgumentException))]
        [TestCase(0, 2, 0, -1, ExpectedException = typeof(ArgumentException))]
        public void GetSubMatrix(int rowStart, int rowLength, int colStart, int colLength)
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix subMatrix = matrix.GetSubMatrix(rowStart, rowLength, colStart, colLength);
                Assert.AreEqual(rowLength, subMatrix.Rows);
                Assert.AreEqual(colLength, subMatrix.Columns);
                for (int i = rowStart, ii = 0; i < rowLength; i++, ii++)
                {
                    for (int j = colStart, jj = 0; j < colLength; j++, jj++)
                    {
                        Assert.AreEqual(matrix[i, j], subMatrix[ii, jj]);
                    }
                }
            }
        }

        [TestCase(0, 2, 0, 2)]
        [TestCase(1, 1, 1, 1)]
        [TestCase(0, 4, 0, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, 2, 0, 4, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(4, 2, 0, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, 2, 4, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(-1, 2, 0, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, 2, -1, 2, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(0, -1, 0, 2, ExpectedException = typeof(ArgumentException))]
        [TestCase(0, 2, 0, -1, ExpectedException = typeof(ArgumentException))]
        public void SetSubMatrix(int rowStart, int rowLength, int colStart, int colLength)
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix subMatrix = matrix.GetSubMatrix(0, 2, 0, 2);
                subMatrix[0, 0] = 1.0;
                subMatrix[0, 1] = 0.0;
                subMatrix[1, 0] = 0.0;
                subMatrix[1, 1] = 2.0;
                matrix.SetSubMatrix(rowStart, rowLength, colStart, colLength, subMatrix);

                for (int i = rowStart, ii = 0; i < rowLength; i++, ii++)
                {
                    for (int j = colStart, jj = 0; j < colLength; j++, jj++)
                    {
                        Assert.AreEqual(matrix[i, j], subMatrix[ii, jj]);
                    }
                }
            }
        }


        [Test]
        public void KroneckerProduct()
        {
            Matrix A = GetMatrix(diag7);
            Matrix B = GetMatrix(diag3);
            Matrix result = A.KroneckerProduct(B);

            for (int i = 0; i < A.Rows; i++)
            {
                for (int j = 0; j < A.Columns; j++)
                {
                    for (int ii = 0; ii < B.Rows; ii++)
                    {
                        for (int jj = 0; jj < B.Columns; jj++)
                        {
                            Assert.AreEqual(result[i * B.Rows + ii, j * B.Columns + jj], A[i, j] * B[ii, jj]);
                        }
                    }
                }
            }
        }


        [Test]
        [ExpectedException(typeof(IndexOutOfRangeException))]
        public void SubtractScalar()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix clone = matrix.Clone();
                matrix.Subtract(2.2);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(clone[i, j] - 2.2, matrix[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(IndexOutOfRangeException))]
        public void SubtractScalarResult()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix result = GetMatrix(matrix.Rows, matrix.Columns);
                matrix.Subtract(2.2, result);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j] - 2.2, result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SubtractScalarResultNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = null;
            matrix.Subtract(2.2, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SubtractScalarResultNotConformable1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = GetMatrix(matrix.Rows + 1, matrix.Columns);
            matrix.Subtract(2.2, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SubtractScalarResultNotConformable2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = GetMatrix(matrix.Rows, matrix.Columns + 1);
            matrix.Subtract(2.2, result);
        }

        [Test]
        public void SubtractMatrix()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix other = matrix.Clone();
                matrix.Subtract(other);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j] - matrix[i, j], matrix[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SubtractMatrixNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = null;
            matrix.Subtract(other);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SubtractMatrixNotConformable1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(matrix.Rows + 1, matrix.Columns);
            matrix.Subtract(other);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SubtractMatrixNotConformable2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(matrix.Rows, matrix.Columns + 1);
            matrix.Subtract(other);
        }

        [Test]
        public void SubtractMatrixResult()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix other = matrix.Clone();
                Matrix result = GetMatrix(matrix.Rows, matrix.Columns);
                matrix.Subtract(other, result);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j] - matrix[i, j], result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SubtractMatrixResultNull1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = null;
            Matrix result = GetMatrix(matrix.Rows, matrix.Columns);
            matrix.Subtract(other, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SubtractMatrixResultNull2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = matrix.Clone();
            Matrix result = null;
            matrix.Subtract(other, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SubtractMatrixResultNotConformable1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(matrix.Rows + 1, matrix.Columns);
            Matrix result = GetMatrix(matrix.Rows, matrix.Columns);
            matrix.Subtract(other, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SubtractMatrixResultNotConformable2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(matrix.Rows, matrix.Columns + 1);
            Matrix result = GetMatrix(matrix.Rows, matrix.Columns);
            matrix.Subtract(other, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SubtractMatrixResultNotConformable3()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(matrix.Rows, matrix.Columns);
            Matrix result = GetMatrix(matrix.Rows + 1, matrix.Columns);
            matrix.Subtract(other, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SubtractMatrixResultNotConformable4()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(matrix.Rows, matrix.Columns);
            Matrix result = GetMatrix(matrix.Rows, matrix.Columns + 1);
            matrix.Subtract(other, result);
        }

        [TestCase(3, 3, 1, 1)]
        [TestCase(3, 3, 1, 2, ExpectedException = typeof(IndexOutOfRangeException))]
        [TestCase(4, 3, 1, 1)]
        [TestCase(4, 3, 1, 2, ExpectedException = typeof(IndexOutOfRangeException))]
        public void Indexer(int row, int column, int i, int j)
        {
            Matrix matrix = GetMatrix(row, column);
            matrix[i, j] = i + j;
            Assert.AreEqual(i + j, matrix[i, j]);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void IndexerGetterArgumentOutOfRange1()
        {
            Matrix matrix = GetMatrix(4, 4);
            double temp = matrix[-1, 0];
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void IndexerGetterArgumentOutOfRange2()
        {
            Matrix matrix = GetMatrix(4, 4);
            double temp = matrix[0, -1];
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void IndexerGetterArgumentOutOfRange3()
        {
            Matrix matrix = GetMatrix(4, 4);
            double temp = matrix[4, 0];
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void IndexerGetterArgumentOutOfRange4()
        {
            Matrix matrix = GetMatrix(4, 4);
            double temp = matrix[0, 4];
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void IndexerSetterArgumentOutOfRange1()
        {
            Matrix matrix = GetMatrix(4, 4);
            matrix[-1, 0] = 1;
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void IndexerSetterArgumentOutOfRange2()
        {
            Matrix matrix = GetMatrix(4, 4);
            matrix[0, -1] = 1;
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void IndexerSetterArgumentOutOfRange3()
        {
            Matrix matrix = GetMatrix(4, 4);
            matrix[4, 0] = 1;
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void IndexerSetterArgumentOutOfRange4()
        {
            Matrix matrix = GetMatrix(4, 4);
            matrix[0, 4] = 1;
        }


        [Test]
        public void ToArray()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                double[,] array = matrix.ToArray();
                Assert.AreEqual(matrix.Rows, array.GetLength(0));
                Assert.AreEqual(matrix.Columns, array.GetLength(1));

                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j], array[i, j]);
                    }
                }
            }
        }

        [Test]
        public void ToColumnWiseArray()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                double[] array = matrix.ToColumnWiseArray();
                Assert.AreEqual(matrix.Rows * matrix.Columns, array.Length);

                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j], array[j * matrix.Rows + i]);
                    }
                }
            }
        }

        [Test]
        public void ToRowWiseArray()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                double[] array = matrix.ToRowWiseArray();
                Assert.AreEqual(matrix.Rows * matrix.Columns, array.Length);

                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j], array[i * matrix.Columns + j]);
                    }
                }
            }
        }

        [Test]
        public void OperatorPlus()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix result = +matrix;
                Assert.IsFalse(ReferenceEquals(matrix, result));
                Assert.AreEqual(matrix.Rows, result.Rows);
                Assert.AreEqual(matrix.Columns, result.Columns);
                for (int i = 0; i < result.Rows; i++)
                {
                    for (int j = 0; j < result.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j], result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatorPlusNull()
        {
            Matrix matrix = null;
            Matrix result = +matrix;
        }

        [Test]
        public void OperatorAddMatrixMatrix()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix other = matrix.Clone();
                Matrix result = matrix + other;
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j] + other[i, j], result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatorAddMatrixNullMatrix()
        {
            Matrix matrix = null;
            Matrix other = GetMatrix(diag3);
            Matrix result = matrix + other;
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatorAddMatrixMatrixNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = null;
            Matrix result = matrix + other;
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void OperatorAddMatrixMatrixNotConformable1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(diag7);
            Matrix result = matrix + other;
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void OperatorAddMatrixMatrixNotConformable2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(diag5);
            Matrix result = matrix + other;
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void OperatorAddMatrixMatrixNotConformable3()
        {
            Matrix matrix = GetMatrix(diag7);
            Matrix other = GetMatrix(diag5);
            Matrix result = matrix + other;
        }

        [Test]
        public void OperatorNegate()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix result = -matrix;
                Assert.IsFalse(ReferenceEquals(matrix, result));
                Assert.AreEqual(matrix.Rows, result.Rows);
                Assert.AreEqual(matrix.Columns, result.Columns);
                for (int i = 0; i < result.Rows; i++)
                {
                    for (int j = 0; j < result.Columns; j++)
                    {
                        Assert.AreEqual(-matrix[i, j], result[i, j]);
                    }
                }
            }
        }


        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatorNegateNull()
        {
            Matrix matrix = null;
            Matrix result = -matrix;
        }

        [Test]
        public void OperatorSubtractMatrixMatrix()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix other = matrix.Clone();
                Matrix result = matrix - other;
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j] - other[i, j], result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatorSubtractMatrixNullMatrix()
        {
            Matrix matrix = null;
            Matrix other = GetMatrix(diag3);
            Matrix result = matrix - other;
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatorSubtractMatrixMatrixNull()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = null;
            Matrix result = matrix - other;
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void OperatorSubtractMatrixMatrixNotConformable1()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(diag7);
            Matrix result = matrix - other;
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void OperatorSubtractMatrixMatrixNotConformable2()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix other = GetMatrix(diag5);
            Matrix result = matrix - other;
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void OperatorSubtractMatrixMatrixNotConformable3()
        {
            Matrix matrix = GetMatrix(diag7);
            Matrix other = GetMatrix(diag5);
            Matrix result = matrix - other;
        }

        [Test]
        public void OperatorDivide()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix result = matrix / .1234;
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j] / .1234, result[i, j], Constants.AcceptableError);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatorDivideNull()
        {
            Matrix matrix = null;
            Matrix result = matrix / .1234;
        }

        [Test]
        public void PointwiseMultiplyResult()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix other = matrix.Clone();
                Matrix result = matrix.Clone();
                matrix.PointwiseMultiply(other, result);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j] * other[i, j], result[i, j]);
                    }
                }

                result = matrix.PointwiseMultiply(other);
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j] * other[i, j], result[i, j]);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void PointwiseMultiplyNull()
        {
            Matrix matrix = GetMatrix(diag7);
            Matrix other = null;
            Matrix result = matrix.Clone();
            matrix.PointwiseMultiply(other, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void PointwiseMultiplyResultNull()
        {
            Matrix matrix = GetMatrix(diag7);
            Matrix other = matrix.Clone();
            matrix.PointwiseMultiply(other, null);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void PointwiseMultiplyNotConformable1()
        {
            Matrix matrix = GetMatrix(diag7);
            Matrix other = GetMatrix(diag7.Length + 1, diag7.Length);
            Matrix result = matrix.Clone();
            matrix.PointwiseMultiply(other, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void PointwiseMultiplyNotConformable2()
        {
            Matrix matrix = GetMatrix(diag7);
            Matrix other = matrix.Clone();
            Matrix result = GetMatrix(diag7.Length + 1, diag7.Length);
            matrix.PointwiseMultiply(other, result);
        }

        [Test]
        public void OperatorMultiplyMatrixScalar()
        {
            foreach (double[] data in matrices)
            {
                Matrix matrix = GetMatrix(data);
                Matrix result = matrix * .1234;
                for (int i = 0; i < matrix.Rows; i++)
                {
                    for (int j = 0; j < matrix.Columns; j++)
                    {
                        Assert.AreEqual(matrix[i, j] * .1234, result[i, j], Constants.AcceptableError);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatorMultiplyMatrixScalarNull()
        {
            Matrix matrix = null;
            Matrix result = matrix * .1234;
        }

        [Test]
        public void OperatorMultiplyScalarMatrix()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix result = .1234 * matrix;
            for (int i = 0; i < matrix.Rows; i++)
            {
                for (int j = 0; j < matrix.Columns; j++)
                {
                    Assert.AreEqual(matrix[i, j] * .1234, result[i, j], Constants.AcceptableError);
                }
            }

            matrix = GetMatrix(diag7);
            result = .1234 * matrix;
            for (int i = 0; i < matrix.Rows; i++)
            {
                for (int j = 0; j < matrix.Columns; j++)
                {
                    Assert.AreEqual(matrix[i, j] * .1234, result[i, j], Constants.AcceptableError);
                }
            }

            matrix = GetMatrix(diag5);
            result = .1234 * matrix;
            for (int i = 0; i < matrix.Rows; i++)
            {
                for (int j = 0; j < matrix.Columns; j++)
                {
                    Assert.AreEqual(matrix[i, j] * .1234, result[i, j], Constants.AcceptableError);
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void OperatorMultiplyScalarMatrixNull()
        {
            Matrix matrix = null;
            Matrix result = .1234 * matrix;
        }

        [Test]
        public void InsertColumn()
        {
            Matrix matrix = GetMatrix(3, 4);
            Vector column = GetVector(matrix.Rows);
            for (int i = 0; i < column.Count; i++)
            {
                column[i] = i;
            }

            for (int k = 0; k < matrix.Columns + 1; k++)
            {
                Matrix result = matrix.InsertColumn(k, column);
                Assert.AreEqual(result.Columns, matrix.Columns + 1);
                for (int col = 0; col < result.Columns; col++)
                {
                    for (int row = 0; row < result.Rows; row++)
                    {
                        if (col == k)
                        {
                            Assert.AreEqual(row, result[row, col]);
                        }
                        else
                        {
                            Assert.AreEqual(0, result[row, col]);
                        }
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void InsertNullColumn()
        {
            Matrix matrix = GetMatrix(3, 3);
            matrix.InsertColumn(0, null);
        }

        [TestCase(-1, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(5, ExpectedException = typeof(ArgumentOutOfRangeException))]
        public void InsertColumnOutOfRange(int row)
        {
            Matrix matrix = GetMatrix(3, 3);
            Vector column = GetVector(3);
            matrix.InsertColumn(row, column);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void InsertNonConformableColumn()
        {
            Matrix matrix = GetMatrix(3, 3);
            Vector column = GetVector(4);
            matrix.InsertColumn(0, column);
        }

        [Test]
        public void InsertRow()
        {
            Matrix matrix = GetMatrix(3, 4);
            Vector row = GetVector(matrix.Columns);
            for (int i = 0; i < row.Count; i++)
            {
                row[i] = i;
            }

            for (int k = 0; k < matrix.Rows + 1; k++)
            {
                Matrix result = matrix.InsertRow(k, row);
                Assert.AreEqual(result.Rows, matrix.Rows + 1);
                for (int r = 0; r < result.Rows; r++)
                {
                    for (int col = 0; r < result.Columns; r++)
                    {
                        if (r == k)
                        {
                            Assert.AreEqual(col, result[r, col]);
                        }
                        else
                        {
                            Assert.AreEqual(0, result[r, col]);
                        }
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void InsertNullRow()
        {
            Matrix matrix = GetMatrix(3, 3);
            matrix.InsertRow(0, null);
        }

        [TestCase(-1, ExpectedException = typeof(ArgumentOutOfRangeException))]
        [TestCase(5, ExpectedException = typeof(ArgumentOutOfRangeException))]
        public void InsertRowOutOfRange(int column)
        {
            Matrix matrix = GetMatrix(3, 3);
            Vector row = GetVector(3);
            matrix.InsertRow(column, row);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void InsertNonConformableRow()
        {
            Matrix matrix = GetMatrix(3, 3);
            Vector row = GetVector(4);
            matrix.InsertRow(0, row);
        }

        [Test]
        public void NormalizeRowsTest1Norm()
        {
            Matrix matrix = GetMatrix(diag3).NormalizeRows(1);
            for (int i = 0; i < matrix.Rows; i++)
            {
                Vector row = matrix.GetRow(i);
                Assert.AreEqual(1.0, row.PNorm(1), 10e-12);
            }
        }

        [Test]
        public void NormalizeRowsTest2Norm()
        {
            Matrix matrix = GetMatrix(diag3).NormalizeRows(2);
            for (int i = 0; i < matrix.Rows; i++)
            {
                Vector row = matrix.GetRow(i);
                Assert.AreEqual(1.0, row.PNorm(2), 10e-12);
            }
        }

        [Test]
        public void NormalizeColumnsTest()
        {
            Matrix matrix = GetMatrix(diag3).NormalizeColumns(2);
            for (int j = 0; j < matrix.Columns; j++)
            {
                Vector col = matrix.GetColumn(j);
                Assert.AreEqual(1.0, col.PNorm(2), 10e-12);
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void NormalizeRowsTestFail()
        {
            Matrix matrix = GetMatrix(3, 3);
            matrix.NormalizeRows(-5);
        }

        [Test]
        [ExpectedException(typeof(ArgumentOutOfRangeException))]
        public void NormalizeColumnsTestFail()
        {
            Matrix matrix = GetMatrix(3, 3);
            matrix.NormalizeColumns(-5);
        }

        [Test]
        public void MatrixMultiplyReturnSelf()
        {
            Matrix matrix = GetMatrix(diag3);
            matrix.Multiply(matrix, matrix);

            Assert.AreEqual(3, matrix.Rows);
            Assert.AreEqual(3, matrix.Columns);
            Assert.AreEqual(-1.1 * -1.1, matrix[0, 0], Constants.AcceptableError);
            Assert.AreEqual(0.0, matrix[0, 1], Constants.AcceptableError);
            Assert.AreEqual(0.0, matrix[0, 2], Constants.AcceptableError);
            Assert.AreEqual(0.0, matrix[1, 0], Constants.AcceptableError);
            Assert.AreEqual(-2.2 * -2.2, matrix[1, 1], Constants.AcceptableError);
            Assert.AreEqual(0.0, matrix[1, 2], Constants.AcceptableError);
            Assert.AreEqual(0.0, matrix[2, 0], Constants.AcceptableError);
            Assert.AreEqual(0.0, matrix[2, 1], Constants.AcceptableError);
            Assert.AreEqual(-3.3 * -3.3, matrix[2, 2], Constants.AcceptableError);
        }

        [Test]
        public void MatrixAddReturnSelf()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix orig = GetMatrix(diag3);
            matrix.Add(matrix, matrix);

            Assert.AreEqual(3, matrix.Rows);
            Assert.AreEqual(3, matrix.Columns);
            Assert.AreEqual(orig[0, 0] + orig[0, 0], matrix[0, 0], Constants.AcceptableError);
            Assert.AreEqual(orig[0, 1] + orig[0, 1], matrix[0, 1], Constants.AcceptableError);
            Assert.AreEqual(orig[0, 2] + orig[0, 2], matrix[0, 2], Constants.AcceptableError);
            Assert.AreEqual(orig[1, 0] + orig[1, 0], matrix[1, 0], Constants.AcceptableError);
            Assert.AreEqual(orig[1, 1] + orig[1, 1], matrix[1, 1], Constants.AcceptableError);
            Assert.AreEqual(orig[1, 2] + orig[1, 2], matrix[1, 2], Constants.AcceptableError);
            Assert.AreEqual(orig[2, 0] + orig[2, 0], matrix[2, 0], Constants.AcceptableError);
            Assert.AreEqual(orig[2, 1] + orig[2, 1], matrix[2, 1], Constants.AcceptableError);
            Assert.AreEqual(orig[2, 2] + orig[2, 2], matrix[2, 2], Constants.AcceptableError);
        }

        [Test]
        public void MatrixSubtractReturnSelf()
        {
            Matrix matrix = GetMatrix(diag3);
            Matrix orig = GetMatrix(diag3);
            matrix.Subtract(matrix, matrix);

            Assert.AreEqual(3, matrix.Rows);
            Assert.AreEqual(3, matrix.Columns);
            Assert.AreEqual(orig[0, 0] - orig[0, 0], matrix[0, 0], Constants.AcceptableError);
            Assert.AreEqual(orig[0, 1] - orig[0, 1], matrix[0, 1], Constants.AcceptableError);
            Assert.AreEqual(orig[0, 2] - orig[0, 2], matrix[0, 2], Constants.AcceptableError);
            Assert.AreEqual(orig[1, 0] - orig[1, 0], matrix[1, 0], Constants.AcceptableError);
            Assert.AreEqual(orig[1, 1] - orig[1, 1], matrix[1, 1], Constants.AcceptableError);
            Assert.AreEqual(orig[1, 2] - orig[1, 2], matrix[1, 2], Constants.AcceptableError);
            Assert.AreEqual(orig[2, 0] - orig[2, 0], matrix[2, 0], Constants.AcceptableError);
            Assert.AreEqual(orig[2, 1] - orig[2, 1], matrix[2, 1], Constants.AcceptableError);
            Assert.AreEqual(orig[2, 2] - orig[2, 2], matrix[2, 2], Constants.AcceptableError);
        }

        [Test]
        public void MatrixMulitplyVectorReturnSelf()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector vector = GetVector(new double[] { 1.1, 2.2, 3.3 });
            matrix.Multiply(vector, vector);
            Assert.AreEqual(-1.1 * 1.1, vector[0], Constants.AcceptableError);
            Assert.AreEqual(-2.2 * 2.2, vector[1], Constants.AcceptableError);
            Assert.AreEqual(-3.3 * 3.3, vector[2], Constants.AcceptableError);
        }

        [Test]
        public void MatrixLeftMultiplyVectorReturnSelf()
        {
            Matrix matrix = GetMatrix(diag3);
            Vector vector = GetVector(new double[] { 1.1, 2.2, 3.3 });
            matrix.LeftMultiply(vector, vector);
            Assert.AreEqual(-1.1 * 1.1, vector[0], Constants.AcceptableError);
            Assert.AreEqual(-2.2 * 2.2, vector[1], Constants.AcceptableError);
            Assert.AreEqual(-3.3 * 3.3, vector[2], Constants.AcceptableError);
        }
    }
}